# Disponivel no pacote de programas como: postrej.py
from numpy import *
from numpy.random import *
from pylab import *

def Likeli(data, limits, nl):
    n = len(data) # Numero de amostras
    data = array(data)
    (ll,ul) = limits #limites do espaco de params.
    step = (ul-ll)/float(nl)
    res = [] #lista de resultados
    sd = std(data) #DP dos dados
    for mu in arange(ll,ul,step):
        res.append(exp(-0.5*sum(((data-mu)/sd)**2)))      
    lik = array(res)/max(array(res)) # Verossimilhanca   
    return lik
    
def amostra(n, data,plotted=0):
    x=uniform(0,1,n) #suporte
    limits = 0,1
    L=Likeli(data, limits, n)
    fx=6*x*(1-x) # priori,beta(2,2)
    s=compress(L[:len(x)]<fx,x) #Rejeicao
    if not plotted:
        p1 = scatter(x,fx)
        p2 = plot(sort(x),L)
        legend([p1,p2],['Priori', 'Verossimilhanca'])
    return s
    
def eficiencia(vector,n, data):
    l = len(vector)
    prob = l/float(n)
    diff = n-l
    n2 = int(diff/prob)
    vec2 = amostra(n2,data,plotted=1)
    s = concatenate((vector,vec2))
    return s,prob
    
def main():
    n=90000
    data = uniform(0,1,3)
    sample=amostra(n, data)
    s,prob = eficiencia(sample,n, data)
    figure(2) #gera histograma
    hist(s, bins=50,normed=1)
    xlabel('x')
    text(0.8,1.2,'Eficiencia:%s'%round(prob,2))
    ylabel('frequencia')
    title('Posterior: n=%s'% n)
    savefig('rejeicao.png',dpi=400)
    show()
    return s   
    
if __name__ == '__main__':
    main() 
